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ABSTRACT 
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Thermodynamics to determine where losses occur. This method has been shown to be superior to 
conventional heat balance analyzes. A primary contributor to the exergy decrease is the production 
of entropy. Entropy is produced in shock waves and boundary layers, as well as in other regions of 
the flowfield. A method of analysis is presented to quantify the entropy that is produced in the bow 


shock wave and the boundary layer of a space launch vehicle for Mach numbers between 1.1 and 10. 
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I. INTRODUCTION 


As commercial interest in space launch vehicles has 
increased over the last twenty years, there has been an 
increased emphasis on improving the efficiency of vehicles 
used for space launch services. This emphasis will continue 
to receive attention as competition forces the providers of 
launch services to decrease costs. 

This thesis is an initial step in the development of an 
exergy methodology tool for the conceptual, preliminary, and 
detailed design of aerospace vehicles. 

In order to determine the efficiency of a articular 
launch system, it is necessary to identify the various 
sources of efficiency loss and the quantity of efficiency 
loss associated with each source within the system. Although 
energy cannot be created or destroyed, it can be transformed 
into unavailable forms. This thesis will present an 
investigation of one method that may be used to examine the 
efficiency loss from two specific areas of a typical launch 
vehicle. This method uses the concept of available energy or 
exergy as the basis for deterrining the decrease in the 
energy state that occurs in the shock wave and the boundary 
layer of a supersonic vehicle. 

Typical performance analyzes of energy systems use only 


the first law of thermodynamics. This method of analysis is 








based on the principle of conservation of energy, and 
suggests that an energy balance accounting for all energy 
flow into or out of the system will provide a measure of the . 


system efficiency. While the first law approach is extremely 





important in the performance analysis of energy systems, it 
provides only the overall performance or the actual gross 
performance of the system or process. To establish the 
magnitude, location, and type of loss within a system it is | 
necessary to use the Second law of Thermodynamics. The | 
second law accounts for the irreversibility of all real 
processes and thus provides a method of calculating the 
losses associated with entropy production that occurs. 

The exergy method of analysis uses both the First and 
Second laws of Thermodynamics to determine the energy that 
is available at various points in a system. While the exergy 


method of analysis has been applied to conventional closed 





cycle thermal systems such as power generation and 
rafrigeration systems, it has not, to the author's 
knowledge, been applied to the open cycle space launch 
vehicle. 

"Exergy is defined as the work that is available ina | 
gas, fluid, or mass as a result of its nonequilibrium 
condition relative to some reference condition."([Ref. 1l:p. 
33] Since exergy is defined as the avaiiable work in a gas, 


fluid, or mass, ict is a special case of defining available 











energy. Exergy is an explicit property, and therefore can be 


calculated at any point in an energy system. 

The exergy method permits the examination of any process 
in relation to the theoretically most efficient manner that 
the process could be carried out within the environment. An 
exergy analysis allows one t% quantify the loss of 
efficiency in a process due to the loss of the quality of 
the energy. 

Exergy is calculated relative to a reference condition 


by the following general equation: (Ref. 1l:p. 34] 


Po ve 
Exergy = (u-u,) - To(s-s,) + ce one =o 7 





cel 
(z- Zo) ae + x (Hem He) No + (1) 
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where the subscript o denotes the reference condition. The 
terms on the right side of the equation are internal energy, 
change in heat energy, flow work, kinetic energy per unit 
mass, potential energy, chemical energy, and radiation 
energy emission respectively. In the analysis of specific 
systems or processes one or more of the above terms may be 
omitted, and only those relevant to the problem being 
considered need be included. [Ref 1:pp.34-36] 

This thesis focuses on the quantity of exergy decrease as 
a result of entropy production across the shock wave and as 
a result of boundary layer formation. As can be seen from 
equation 1, an increase in entropy results in a decrease in 


3 


exergy. The Second law of Thermodynamics tells us that all 


real process are irreversible and that entropy is increased 


in an irreversible process. 














II. EXERGY DECREASE ACROSS THE SHOCK WAVE 


Supersonic flow over a blunt body will result in the bow 
shock wave being detached from the body as shown in Figure 
1. The shape of the detached shock wave, its detachment 
distance, and the flowfield between the body and the shock 


are determined by the shape of the body and the free stream 


Mach number. [Ref.2:p.128-129] 


Uniform free stream 





Figure 1: Supersonic Blunt Body in Flowfield. 








Since the velocity of the body is supersonic, the flow ahead 
of the body is unaware of the presence of the body until the 
shock wave is encountered. The shock wave directly in front E 
of the body is approximately normal to the free stream and 
consequently the flow behind the shock wave at this point is 
subsonic. The subsonic flow region is defined by the shock 
wave, the body, and the sonic line as shown in Figure 1. 
Because the flowfield between a blunt body and its 
detached shock wave contains both subsonic and supersonic 
flow, an exact determination of the flowfield dynamics is 
extremely complicated to formulate. The steady state 
flowfield in the subsonic region is described by elliptic 
equations ard in the supersonic region by hyperbolic 
equations. The sonic line divides these two regions. It is 
the problem of predicting the flowfield in the transition 
region between the subsonic and supersonic flow that makes 
the flowfield calculations for a supersonic blunt body so 


adifficult. 





The flowfield behind a detached shock wave cannot be 
precisely determined without the application of 
sophisticated numerical methods. In keeping with the 
development of a design tool, a simplified method for 
determining the approximate shape of the detached shock wave 
that was developed by Moeckel [Ref. 3] is used. This method 


was subsequently examined by Love [Ref. 4] and with some 








modification was found to be valid for blunt bodies 
traveling at supersonic and hypersonic speeds. 

Love introduced two major modifications to Moeckel’s 
method. The first of these was in the method used to 
calculate the shock detachment distance from the body. 
Moeckel used a continuity relation to calculate detachment 
distance and shock shape. Love calculates shock detachment 


distance by the following equation: 


x! 


3, = 9-5C cot dace (2) 


where x’/d’ is the nondimensional detachment distance, C is 
a factor that depends on type of body (spherical, 
cylindrical, etc.), and 6ga¢ is the required deflection 
angle for shock wave detachment, from a conical body, for 
the given Mach number (see Figure 2). 

The second modification made by Love was in the 
calculation of the angle eta (7) that the sonic line makes 
with the freestream. Love conducted wind tunnel tests and 
compared his experimental results with his theoretical 
calculations [Ref.4:pp.13-15]. For this analysis, the angle 
” corresponding to 1.1sMs10 was obtained by developing the 
followirg analytical expression for the values of y vs Mach 


number as set forth by Love [(Ref.4:p.45]: 


= 3.9901+56.9298* log(M) (3) 
57.29578 











Figure 2: Representation of Flow With Detached Shock 
and Notation Used in Analysis. 


A. BODY SHAPE 

The shape of the deflecting body is one of the factors 
determining the shock wave shape. For this study, the body 
was patterned somewhat loosely after the Titan III-A launch 
vehicle. This body shape is essentially a blunt cone and 
cylinder arrangement. The dimensions of the launch vehicle 
used in the analysis are shown in Table 1. 


TABLE 1: LAUNCH VEHICLE DIMENSIONS 


10 Ft. 
1.5 Ft. 












Nose Radius 














B. MACH NUMBER 

Besides the shape of the body, the other factor in 
determining the shape of the shock wave is the freestream 
Mach number. This study determined the shock shape based on 
1.1sMs10. Although this is far short of the Mach number 
required to reach orbit, the purpose of the study was to 
demonstrate the validity of the method, which can be 


accomplished in the listed Mach range. 


C. ASSUMPTIONS 
The method of approximating the shock shape and location 
developed by Moeckel is based on two assumptions: 
(1) The form of the shock between its foremost point 
and its sonic point is adequately represented by 
a hyperbola asymptotic to the free stream Mach 
lines; and 
(2) the sonic line between the shock and body is 
straight and inclined at an angle that depends 
only on the freestream Mach number. [Ref.3:p.1] 
The first assumption states that the shock shape between its 
vertex, or foremost point, and its sonic point may be 
represented by a hyperbola. For this study, the assumption 
was made that this approximation for the shock shape is 
valid past the sonic point of the shock wave. For 
calculation purposes, the shock wave was taken from the 
vertex of the wave to the point where the angle between the 
shock wave and the asymptote that defines the hyperbola 


differ by no more than one per cent. The equation for the 


hyperbola used in these assumptions is 








“Xo (4) 





where x is the longitudinal axis of the body, y is the axis 
perpendicular to the x axis, 8 is the cotangent of the Mach 
angle, and x, is the distance from the vertex of the shock 
wave to the intersection of its asymptotes. [Ref.3:p.1] 

The second assumption made by Moeckel concerning the sonic 
line was modified by Love as previously discussed and was 
used in this analyzes. 

Air is the fluid used for all of the analysis presented 
in this thesis. However, the analysis method is equally 
valid for any gas, and would only require using the 
appropriate value for the ratio of specific heats (y) and 
the new specific gas constant. It has been assumed that the 
fluid behaves as a perfect gas in all of the calculations. 
This assumption is valid for velocities up to approximately 
M=3 to M=5S [Ref 5:pp.17-20]. Above these velocities 
dissociation of the fluid begins to occur due to atmospheric 
heating. If the gas continues to be heated, ionization will 
occur. As the degree of dissociation and ionization 
increases, the perfect gas assumption becomes less valid. 

The calculation of total entropy across the shock wave 
requires a determination of the mass flow rate of air 
through the shock wave. The air mass flow rate calculation 
is based on the density, cross sectional area being studied, 
and the velocity of the air. Since the density of air in the 


10 


atmosphere varies with height above sea level, an arbitrary 
height of 30,000 meters was assumed for all M values. The 
density and temperature of air at this altitude was taken 
from the International Civil Aviation Organization (ICAO) 
standard atmosphere tables [Ref 6]. The calculations may be 
used for any altitude by inserting the appropriate density 


and temperature. 


D. EXERGY CALCULATION 

From equation (1), and neglecting all terms on the right 
except the entropy term, the total exergy decrease through 
the shock wave is proportional to the increase in entropy 
across the shock. The entropy change across a shock is given 


by [Ref. 7:p. 10) 


2.3 2 22303 
27M. 6-(y-1 1) ™ 6 
As ee yM, sin (y yy Seal (y+1)M, sin 


ey yrs (y-1) M2 sin? 6 +2 (5) 
where c, is the specific heat at constant volume for air, y 
is the ratio of specific heats, M is the freestream Mach 
number, and @ is the local shock wave angle with reference 
to the freestream direction. By using the method of Moeckel, 
modified by Love, to determine the approximate shock shape, 
the entropy increase across the total shock was calculated 
by equation 5 and numerical integration using cumulative 
applications of Simpson's Rule. The number of integrations 


varied with the shock strength and shape running from a low 


of 50 for Mach=1.1 to a high of 2513 for Mach=10. 


a 











The entropy calculations were performed by the Fortran 
77 program contained in Appendix A. For the assumed body 
shape, this program calculates the approximate shock wave 
shape and the total rate of entropy change across the shock 
for a designated control volume and designated Mach numbers 
between 1.1 and 10. The calculations are carried to three 
decimal places with the knowledge that assumptions and 
round-off error may have affected the accuracy. The results 
are summarized in Table 2. 


TABLE 2: ENTROPY CHANGE ACROSS SHOCK 


Mach Number Rate of Entropy Change 
Across Shock (Kd/sec} 





; 

50.803 

301.469 

a 836.941 

1662.150 
. 











As can be seen from equation 5, the entropy rise across 
the shock increases as the shock wave angle (@), and the 
freestream Mach number (M,) increase. For a detached shock 
wave, the shock wave angle is greatest immediately ahead of 
the body where it is approximately normal to the freestream 


and decreases as it curves around the body until it 
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eventually asymptotes the freestream. This indicates that 

for a given body shape and Mach number, the greatest exergy 
loss occurs in the vicinity of the vertex of the shock wave. 
This is shown graphically in Figure 3. The vertical axis is 
shown as entropy change over the specific gas constant (R’), 
and the horizontal axis is shown as the angle phi (¢), which 


is 90 - @. 
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Figure 3: Entropy Rise Across Shock Wave. 


E. VERIFICATION OF THE COMPUTER MODEL 

Comparisons between theoretical shock location using the 
program of Appendix A and actual shock location as 
determined by Love for Mach numbers of 1.94 and 6.8 fora 
combination sphere and cylinder body shape are shown in 
Figures 4 and 5 respectively. These comparisons show good 


13 








agreement between calculated shock shape and shock location. 
The data points corresponding to experimental results were 
taken from charts [Ref. 4:pp 46,51] and therefore, precise 
quantitative evaluations regarding goodness of fit are not 
possible. However, the experimental data appear to be within 


five percent of the predicted values of the shock profile. 


2m} /q? 
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Ss i ae 
ee 

sree eneaae 
7 a 


(MXN gd 





Figure 4: Comparison of Experimental and Predicted 
Shock Shape. Mach 1.94 





Figure 5: Comparison of Experimental and Predicted 
Shock Shape. Mach 6.8 
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Another reality check for the computer model was made by 
calculating the mass flow of air through the control volume 
defined by the shock wave at 1.1sMs10 at an altitude of 
30,000 meters. Since the vehicle is supersonic, there is no 
deflection of the flow around the body prior to passing 
through the shock. The validity of the shock shape generated 
by the program is good as shown in Figures 4 and 5 above. 
The shock wave may be thought of as a blunt cone. If the 
cone is terminated at any point, the distance from the 
centerline to the shock is the radius of the base. The mass 
flow through the area of the base can be easily calculated. 
This was done and the result was compared with the mass flow 
calculated by the program using the same Simpson's rule 
integration technique that was used to calculate the entropy 
rise across the shock. The results of these two calculations 


differed by less than three per cent for all Mach numbers. 
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IIIT. EXERGY DECREASE IN THE BOUNDARY LAYER 


Tf a control volume is taken around a cylindrical body 
moving through a gas, there wili be a change in the momentum 
of the gas taken upstream of the body from that taken 
downstream of the body. This difference in momentum is a 
function of the drag on the body, and is given by [Ref. 


8:pp. 23-24] 
D= f pu(u-u)da (6) 
A 


where the integral is taken over the end surface of the 
cylinder at a large distance downstream, U is the velocity 
upstream, and u is the velocity at an arbitrary point in the 


wake. Letting 


tie dif fat 0 
as Lge 5302) (7) 


According to Crocco’s theorem ‘Tef. 9] 

- (Hays Maz) = T Bay+ Baz) (8) 
ly Zz u oy Zz 

Substituting equations 7 and 8 into equation 6, we get 


= ° TOs ds 
D= [ puaal. 7 Say ae (9) 


Replacing pudA (the mass flow through the cross-sectional 
area GA) by dm and T/u by T,/U, which is justified if second 


order terms are neglected [Ref. 8:p.24], in equation 9 gives 
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Ds = [ dm(s-s.) (10) 
4 


Taking the integral over the area of the wake, equation 10 


can be written as 


DU = T= (11) 
sec 


This is an application of the theory of Oswatitsch [Ref. 
10]. The freestream temperature and velocity are known or 
can be calculated using the ICAO tables [Ref. 6] for a given 
altitude and Mach number. If the drag can be calculated, 
then the entropy increase can. be calculated using equation 
11. To calculate the profile drag of the body, the flow 
parameters at the boundary layer edge and the boundary layer 
shape are needed. 

The theory that the flow surrounding an immersed body 
could be divided into two parts was first proposed by L. 
Prandtl in 1904. Prandtl stated that the viscosity in the 
major part of the flow field has a negligible effect on its 
motion and consequently in most cases could be ignored. The 
other part of the flow field consisted of a narrow region 
next to the body where the viscous forces are not 
negligible. [Ref.11] 

This narrow region near the body is known as the 
boundary layer. In the boundary layer the velocity of the 
fluid increases rapidly from zero at the surface to a value 


which then changes slowly with increasing distance from the 
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body. At an infinite distance from the body, the boundary 
layer and the freestream merge smoothly together. For 
practical applications, the boundary layer conditions 
approach the freestream condition very rapidly within a 
relatively small distance from the surface as shown in 


Figure 5. 





Figure 6: Boundary Layer on a Flat Plate [Ref. 10:p. 25] 


If the freestream velocity is taken as U, then the edge of 
the boundary layer may be designated as that distance from 
the body where the velocity of the stream is .995U or u,. 
Boundary layers may be either laminar or turbulent. In 
laminar flow, the flow within the boundary layer is 
characterized by smooth streamlines approximately parallel 
to the surface and there is little mixing. In the turbulent 
boundary layer the flow experiences irregular fluctuations 
in velocity magnitude and direction superimposed on the mean 
flow, which remains roughly parallel to the surface 
{Ref. 12:p.11]. 
The change from laminar to turbulent flow is referred to 


as transition. Transition may be caused by a number of 
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factors including disturbances that occur in the freestream, 





surface roughness, changes in the pressure distribution on 


the surface of the body, or an increase in Reynolds number. 

























Reynolds number is a ratio of the inertia to viscous forces 
in a fluid and is given by the following equation: 


Re sear (12) 


where U is the velocity, p is the density, 1 is the body 
dimension, and p is the dynamic coefficient of viscosity of 
the fluid. 

The presence of shear forces in the boundary layer due 
to friction and the viscous nature of the fluid result in 
heat being generated. The heat is transferred to the body 
and to the fluid. The viscous dissipation as heat in the 
boundary layer results in an increase in entropy and 
consequently a decrease in exergy. At hypersonic velocities, 


che intense viscous dissipation may result in temperatures 





high enough to cause dissociation and ionization of the 
surrounding air [Ref. 5:p. 228]. This phenomenon and its 
effects will not be discussed in this paper. However, the 
error introduced by ignoring these effects can be 
significant. Interpolating from Anderson [Ref. 5:p. 511], 
the error, at M=10 and at 30,000 meters, in the ratio T,/T, 


across a normal shock is approximately 24 per cent. 
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A. BOUNDARY LAYER CALCULATIONS 
An exact determination of the flow conditions in and 
surrounding the boundary layer for a body of arbitrary shape . 
is possible by solving the three dimensional Navier-Stokes 
equations. These equations do not have a closed form 
solution, but require elaborate computer methods to be 
solved. However, if we specify that the boundary layer 
thickness (6) is much less than the length of the body, and 
that the Reynolds number is large (>10°), the Navier-Stokes 


equations can be reduced to the boundary layer equations: 





Continuity: G(pu) , O(pv) LG (13) 
Ox oy 

; du,,,,du . _ We, 2, du 4 

x Momentum: Uae Pay ae + By (H sy) ( ) 


uN 
Oo 


y Momentum: 2p (15) 


; oh Oh_ 9 ,,9T ap. 
Energy: pu Pay Tae as ag 





+H Say? (16) 


where x is the direction parallel to the surface, y is the 
direction normal to the surface, u is the velocity component 
in the x direction, v is the velocity component in the y 
direction, p is the density, p, is the pressure at the outer 
edge of the boundary layer, nw is the viscosity ccefficient 
of the fluid, h is the enthalpy, and k is the thermal 


conductivity of the fluid. 
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Like the Navier-Stokes equations, the boundary layer 
equations are nonlinear. They are simpler and easier to 
solve, but still require computer methods to obtain exact 
solvtions in most cases. Solution methods for the boundary 
layer equations are presented by Young [Ref. 12] and 
Schlichting [Ref. 14]. 

1. Approximate Methods 

The purpose of this thesis was not to determine an 
exact solution to the boundary layer equations, but to show 
that an exergy analysis can be conducted on the boundary 
layer and to provide the foundation for an exergy based 
design tool and or methodology. This concept can be 
demonstrated by making several approximations that will 
' provide rough estimates of the required variables without 
the extensive numerical and computer calculations required 
to obtain a more exact solution of the boundary layer 
equations. The necessary calculations are developed in the 
following paragraphs. 

a. Stagnation Point Calculation 

As stated in the preceding chapter, the flow 
field surrounding a blunt body traveling at supersonic 


velocity will consist of both subsonic and supersonic flow 


divided by the sonic line (see Figure 1). For a body at zero 


angle of attack, the flow precisely along the body axis will 


pass through the shock wave, which at that point is normal 


or very nearly normal to the freestream direction, slow, and 
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stop at the stagnation point. On a blunt body, the boundary 
layer has a finite thickness at the stagnation point. The 
boundary layer thickness at the stagnation point of a two 
dimensional body is given by 


8(0) = ote (17) 
e 


where 6(0) is the boundary layer thickness at the stagnation 
point, and V is the magnitude of the velocity gradient at 
the stagnation point. [Ref. 15:pp. 208-209] 

The velocity gradient is determined by the 
method of Anderson [Ref. S:pp. 255-256]. It can be shown 
using Euler’s equation applied at the boundary layer that, 


du, 1 ap, 


= TE (18) 











where the subscript e indicates conditions at the edge of 
the boundary layer, u is the velocity component in the x 
direction which is parallel to the body, and p is the 
pressure. Assuming a modified newtonian pressure 
distribution over the surface, the pressure coefficient is 


given by [Ref. S:p. 53] 


C. = 


> = Cp sin’@, (19) 


where 6, ig the angle between the tangent to the surface and 
the freestream direction, and Comax is given in terms of y 


and the freestream Mach number (M,) as 
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(y+1)2M2 Y. 4-y+3yM? 
7 Diet en ease i pa alae ike a pee (20) 


pee 
4yM? -2(y-1) vo 
Defining ¢ as the angle between the normal to the surface 


and the freestream, then equation 19 can be written as 


Cy = Cy .cOs* > (21) 
Defining Cp as 
C, Pew Pi (22) 
Y 


where the subscript 1 denotes the freestream condition, then 


equation 21 can be written as: 


P, “Py; 


3 = CoC OS" O (23) 
1 


where q, is the dynamic pressure of the freestream. Equation 


23 can be written as 
Pe = Pang Ft COS? O +P, (24) 
Differentiating equation 24 with respect to x, we obtain 


dp, 
dx 





Sa ng oO 
= 2c, .4, cos sing — (25) 


Substituting equation 25 into equation 18, we have 





du, 2c, q, < a 
= —Poax 2} ab (26) 
ax 3.u, cososing dx 


Consider the stagnation region, as shown in Figure 7. 
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Figure 7: Stagnation Region of a Blunt Body [(Ref.5:p. 
251). 


In this region, let Ax be a small increment of surface 
distance above the stagnation point, corresponding to small 
changes in ¢ of Ad. Since the flow velocity is so low in the 


stagnation region, we can assume [Ref. 5:p. 252] 


du, 
ax ) pas (27) 





u, = ( 


In the stagnation region the following approximations can be 


made if ¢ is small: 


cos = 1 (28a) 
; he bAn AX 
sing =o=Ag=—- (28b) 
ao 2 1 (28¢) ; 
ax R 


where R ig the radius of the nose of the body at the 


stagnation point. 
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With the approximation of equation 28a and 


substituting into equation 22, equation 21 becomes 


P,-P 
Ch = Coax =, (Soe eek (29) 
ay 
which is rewritten as 
= 1 
Gq = Go (Pe-P;) (30) 
Prax 


Substituting equations 27-30 into equation 26 and 


rearranging terms gives 


du 


e . 1:| 2(P.-P,) (31) 
ax R P, 


This value for the velocity gradient is then substituted 





into the value for V in equation 17. 

The velocity, temperature, density, Mach number, 
and pressure immediately behind the shock wave can be 
determined from the normal and oblique shock relations 
contained in Anderson [Ref. 2:pp.67-68,106]. The density 
behind the shock along the stagnation streamline is used in 
equation 17 for the value of p. 

Equation 17 gives the boundary layer thickness 
at the stagnation point of a two dimensional body. The 
boundary layer thickness for an axisymmetric body will be 
thinner than for the two dimensional case [Ref. 5:p. 254]. 
While the difference between the axisymmetric and two 


dimensional cases is noted, equation 17 does provide an 
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approximation of the boundary layer thickness for the 
axisymmetric case at the stagnation point. 

b. Boundary Layer Edge Conditions 

The pressure distribution at the surface of the 

body is approximated by the modified newtonian method (Ref. 
5:pp. 53-56]. Although the results provided by this method 
are not very accurate at low supersonic Mach numbers, the 
results improve quickly as Mach number increases [Ref.5:pp. 
45-56]. The pressure coefficient at the surface of the body 
using the modified newtonian method is given by equation 19. 
Recalling from equation 22 the definition for Cp, the 
pressure at the boundary layer edge over the nose of the 


body is approximated by 
Pe = Cp 900876 +p, (32) 


Figure 8 shows the gurface pressure distribution over the 
nose of the hypothetical vehicle used in this analysis for 
selected Mach numbers. The nose radius (R) in Figure 8 is 
1.5 ft (.4572 m) and the pressure distribution is shown from 
the stagnation point through an angle of 68.2 degrees, where 
the rounded nose intersects the cone of the body. 

The total temperature behind the shock wave 
remains constant everywhere behind the shock. Also, we 
assume that the total pressure is constant along the 
streamline that enters through the nearly normal portion of 


the shock and flows along the boundary layer edge. These 
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values can be calculated and used as a reference to obtain 


local values of temperature and pressure. 
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Figure 8: Pressure Distribution Over a Blunt Nose 
Cylindrical Cone Body. 


The total pressure in the freestream is given by 





at oe 
Po, = P[(1+ yet) y-1] (33) 


where P,, is the stagnation pressure in the freestream, P is 
the static pressure of the freestream. For these analyzes 
all freestream atmospheric conditions were taken from 
Reference [6] at an altitude of 30,000 meters. The ratio of 


stagnation pressures across the shock is 
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rar ae (34) 


where R’ is the specific gas constant, and s,-s, is given by 








2 
2+ (y-1)M 
5-8 Senin ae ay) 
yr (y +1) M2) (35) 
is + 27 (mM? - 
Rin(1+ 2% (My -1)] 


where Cc, is the specific heat at constant pressure of the 
gas. The stagnation pressure behind the shock is then 
obtained by taking the product of equations 33 and 34 


oe ae 


Po, By 
1 


) Po, (36) 


The stagnation temperature in the freestream is 


calculated by 


1 


T,=T(1+ 1 M;) (37) 





where T is the freestream temperature. The total temperature 


across a normal shock is constant, therefore 


= To (38) 


If we substitute Ty), for Tp; in equation 37, and with a 
known Mach number, the local temperature can be determined. 
The sonic point on the body nose is located with 


the same method used in determining the shock shape, and is 


described in Appendix A. Since the Mach number at the sonic 





line is one, the local pressure at the sonic point can be 
obtained by substituting the value of Py),, from equation 36, 
for Py, in equation 33. Once the sonic point is located, the 
angle between the body centerline and the sonic point can be 
easily calculated. This angle will be designated the sonic 
point angle (SPA) for calculation purposes. The total angle 
that the flow makes as it rounds the nose is 68.2 degrees. 
Subtracting the SPA from 68.2 gives the remaining angle, 
designated A,, that the flow must turn through in traveling 
from the sonic line to the cone portion of the body. 
Treating angle A, as a Prandtl-Meyer expansion angle, the 
pressure, Mach number, and temperature of the flow 
downstream of the expansion can be calculated. 

Before proceeding further, it is necessary to 
introduce the Prandtl-Meyer function [Ref. 2:p. 134], which 


is given the symbol p 





vy(M) = ES tan*| 22-1) -tan7+¥M?-1 (39) 
y-1 yr1 


As can be seen from equation 39, the Prandtl-Meyer function 
is a function of the Mach number and y. Let v(M,) be the 


Prandtl-Meyer function before the expansion, and »(M,) be | 





the Prandtl-Meyer function after the expansion. Then the 


angle A, that the flow turns through is 
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A; = vy (M>) ~v(M,) (40) 


At the sonic line the Mach number is one, and 
therefore »(M,) in equation 40 is zero and A,:»(M,). Putting 
the known value of A, into equation 39 for »(M), the Mach 
number after the expansion can be determined. This process 
is repeated at the point wh?re the cone joins the body. In 
this fashion the flow Mach number along the length of the 
cylindrical portion of the vehicle can be found. With the 
Mach number known and the previously calculated values of 
total pressure and total temperature, the local pressure and 
temperature along the cylindrical portion of the body can be 
determined using equations 32 and 37. These values will be 
taken as the conditions at the edge of the boundary layer. 

c. Boundary Layer Thickness 

The boundary layer is formed as a result of the 
no slip condition that exists at the wall of a body immersed 
in a flow (see Figure 6). As the flow progresses along the 
length of the body, the thickness (6) of the layer of 
retarded flow increases as more and more fluid becomes 
affected (Ref. 14:p. 25]. The thickness of the boundary 
layer on a flat plate at zero incidence can be calculated 
using the method developed by Blasius [Ref. 16]. For laminar 


flow the boundary layer thickness was found to be 


(41) 
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where x is the distance along the plate, and Re, is the 
Reynolds number (see equation 12) as a function of distance 
along the plate. For a turbulent boundary layer the 


thickness is given by [Ref. l2:p. 13] 


0.37x 


6= 
(Re) 1/5 (42) 


As previously stated, the boundary layer on an 
axisymmetric body is thinner than on a flat plate. However, 
as a first approximation, the boundary layer over a flat 


plate with a length equal to the distance traveled by the 





flow from the stagnation point to the tail end of the body 
will be used to approximate the boundary layer thickness. 
This approximation is reasonable when 6<<R as it is in this 
case. 

Since there is a dramatic difference in the 
growth of laminar vs turbulent boundary layers, it is 
necessary to establish the point where transition occurs. In 
actuality transition may occur over a finite distance making 
prediction of a point of transition extremely difficult. 
Research to determine exactly how transition occurs and 
methods to predict its occurrence are ongoing. A formula 
developed by Anderson [Ref. 5:p 280] to predict the point of 


transition will be used in this analysis. 
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logy (Rez) = 6.421exp[1.209x10-4m2-&4?) (43) 


where Re, is the Reynolds number at which transition occurs, 
and M, is the Mach number at the boundary layer edge. It 
should be recognized that the data used to develop this 
equation was unknown to the author and that equation 42 was 


being used as an approximation. 


B. BASE PRESSURE DETERMINATION 

If the assumption is made that the launch vehicle is a 
projectile with no exhaust plume, the pressure behind the 
body in the wake is assumed to be the freestream static 
pressure. This assumption is supported by Anderson [Ref. 
5:pp 117-138] who shows that for a cylindrical body with 
fineness ratio greater than six, the pressure asymptotically 
approaches freestream pressure. The assumption of no exhaust 
plume allows calculation of the entropy change as a result 
of the boundary layer on the body, but ignores the effects 
of boundary layer/plume interaction. 

With the pressure in the wake being approximately equal 
to static pressure and the stagnation pressure behind the 
shock, from equation 36, the Mach number in the wake can be 
determined using equation 33. The difference between the 
Mach number in the wake and in the flowfield can be used to 
calculate the difference between the velocities (U-u) needed 


to solve equation 6. These values are summarized in Table 3. 
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C. WAKE PROFILE 

The wake boundary is established from the body plus the 
boundary layer thickness at the base of the body. Traveling 
downstream for several body diameters and assuming the same 
rate of growth of the boundary layer, the distance from the 
centerline to the point where the wake velocity approaches 
the freestream velocity can be determined. This is shown in 
Figure 9. For calculation purposes, a distance of five 
diameters downstream of the base was arbitrarily chosen. 


TABLE 3: COMPARISON OF FREESTREAM AND WAKE MACH NUMBER 
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Figure 9: Illustration of Wake Profile. 
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The variation of u in the wake from the centerline to 
the point where u*U is given by 
U, 
u=U,+ x(1~—) JU, (44) 
U; 
where U, is the minimum wake velocity, U, is the freestream 
velocity, and x is a scaling variable where Osxsl. The 
distance from the centerline to the point where u-U is given 
by 
y = (R, + 8.) 8 (45) 
where R, is the radius of the base of the body and 6, is the 


thickness of the boundary layer extended to the point in the 


wake where the profile is being calculated. 


D. DRAG DETERMINATION 
The drag of the body is calculated using equation 6 


which is repeated here: 
D = f pu(u-u)da (46) 
A 


Since the velocity in the wake varies along the centerline, 
designated the x axis, and with the distance from the 
centerline, designated the y axis, equation 46 becomes the 


following double integral: 
D= ["[(2xzeu(u-u) dr (47) 
Qo 40 


where u and y are given by equations 44, and 45 
respectively. The density in equation 47 is the density at 
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the location in the wake where the profile is determined. 
For these calculations this value was approximated using the 
perfect gas equation. The pressure is taken as freestream 
and the temperature is that calculated for the boundary 
layer temperature just forward of the base of the body 


ignoring plume effects. 


E. ENTROPY CALCULATION 
With the drag determined, the entropy flux increase as a 
result of the boundary layer can now be calculated using 


equation 11 which is repeated here 


DU=T, = (48) 





where U and T, are the freestream velocity and temperature 
respectively, and D is the drag. 

For the body shape used in this analysis and with the 
approximations used to calculate the required flowfield 
Parameters, the entropy increase is presented in Table 4. 

As can be seen from Table 4, the entropy increase in the 
boundary layer starts to decrease above M=7. This occurs due 
to the boundary layer becoming thinner at higher Ma 1 
numbers. The thinner boundary layer is a result of the 
transition point moving aft on the body. Initially the 
transition point moves forward toward the nose of the body. 
Between 6sMs7, the transition point stops moving toward the 
nose and starts moving back toward the tail. This is 
occurring as a result of the laminar to turbulent 
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TABLE 4: ENTROPY INCREASE DUE TO THE BOUNDARY LAYER 


Mach number Rate of Entropy Increase 
{KJ /sec) 





transition point as calculated by equation 43. The 
transition Reynolds number increases exponentially with 
increasing Mach number. For a given transition Reynolds 
number, the point of transition is obtained by substituting 
Re, for Re and x for 1 in equation 12, and rearranging 


Re rH 
pV 


(49) 





As Re, increases exponentially, the numerator in equation 49 
grows faster than the denominator, as shown in Figure 10. 
The result of this is that the transition point moves toward 
the tail of the body and the boundary layer remains laminar 
longer. Since entropy is produced at a higher rate in 
turbulent boundary layers, the net result is that less 


entropy production occurs. 
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Figure 10: Variation of Transition Point Parameters With 
Mach Number. 


This phenomenon is a result of using equation 43 to 
predict the transition point and may not be an accurate 
representation of actual flow. As previously stated, the 
data that equation 43 was developed from was unknown to the 
author, and the conditions of this study may be outside the 


range for which equation 43 is valid. 
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Iv. CONCLUSIONS 


Recall from the exergy equation (equation 1) that an 
increase in entropy results in a corresponding decrease in 
exergy. A primary advantage of exergy analysis is the 
identification of those processes or areas that have the 
highest exergy loss and hence the greatest potential for 
improved efficiency. A method to quantify the exergy 
decrease resulting from entropy production across the shock 
wave and in the boundary layer, ignoring high temperature 
effects, has been presented. 

It has been demonstrated that the exergy decrease due to 
the production of entropy across the shock wave and in the 
boundary layer can be calculated for a typical space launch 
vehicle. This information can be used to highlight those 
areas which have the greatest potential for improved 
efficiency. 

For the conditions of this study and within the accuracy 
of the approximations, it is clearly evident that entropy 
production +s much greater in the boundary layer than across 
the shock wave for a given Mach number. This indicates that, 
all other factors being equal, design efforts should 


concentrate on reducing the entropy production in the 


boundary layer. 











V. RECOMMENDATIONS FOR FURTHER RESEARCH 


This thesis applies the method of exergy analysis to two 
areas of interest for space launch vehicles. In doing so, it 
provides an initial step in the development of an exergy 
methodology tool for design purposes. 

Following the work of Czysz [Ref. 21], application of 
this methodology to launch vehicle propulsion systems would 
be beneficial. 

Alternately, variations in launch vehicle design could 
be examined with the goal of finding the optimum design to 
minimize entropy production across the shock wave and within 
the boundary layer. 

It would also be useful to determine the magnitude of 
the error for the approximations used in this study. This 
could be accomplished by comparing the results with those 
obtained from exact flowfield calculations using CFD 
methods. Also the effects of the boundary layer/plume 


interactions on entropy production could be investigated. 
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APPENDIX A 
THE SHOCK WAVE MODEL 


This program computes the entropy increase across the 


shock wave of a cylindrical space launch vehicle with a 


blunt cone nose. The program is written for Mach numbers of 


1.1 to 10, and for atmospheric conditions at an altitude of 


30,000 meters. This program is written in Fortran 77 code. 


Comment statements are preceded by a "C", 


qaaagagNnaaNnaAaaNnanaNnaaaAaNnANnanaANAANAANanaa a 


Variable declarations 


M = Freestream Mach number 

BETA = Cotangent of the Mach number 

THETA = The shock wave angle for sonic flow behind the 
shock 

DELTA = The maximum stream deflection angle that can 
occur without shock separation 

NU = Angle between the sonic line and a line normal to 
the freestream direction 

GAMMA = Ration of specific heats (1.4 for air) 

SIGMA = Angle of shock wave at the center of mass 
streamline 

A = Area ratio for isentropic flow 

P = Stagnation pressure ratio 

RATIO = Ratio of the distance from the x-axis to the 
shock sonic point to the distance from the x-axis 
to the body sonic point 

YSB = Distance from the centerline of the body to the 
sgonic point on the body 

XO = Distance from foremost point of detached shock to 
intercept of its asymptote on the x-axis 

Y = Coordinate perpendicular to freestream flow 

X = Coordinate in the direction of the freestream 

I,J,%2 = Increment counters 

PI = Radian equivalent of 180 degrees 

DIFF = The difference between the shock wave angle and 
the asymptote angle of the hyperbola that is 
defined by the shock wave 

CANGLE «= The constant angle between the asymptote of 
the hyperbola and the freestream direction 
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SANGLE = The angle between the freestream direction and 
the tangent to the shock wave 

ANGLED = The game angle as "SANGLE" but in degrees 
instead of radians 

TOTALS = The total change in entropy across the shock 
wave 


START PROGRAM 

INTEGER I, J, Z 

PARAMETER (2Z=10) 

REAL M(Z), BETA(Z), THETA(Z), DELTA(2), 


+ NU(Z), PI, YSB(Z), SIGMA(Z), A(Z), P(Z), 
+ RATIO(Z), X0(Z), X(7500), Y¥(7500), DIFF, B(Z), 
+ CANGLE(Z), SANGLE(7500), ANGLED(7500), TOTALS 
PI=3.14159 

M(1)=1.1 

M(2) =2 


BETA (1) =SQRT (M(1) **2.0-1.0) 
BETA (2) =SQRT (M(2) **2.0-1.0) 
DO 05 I=3,2 
M(I) =M(I-1)+1.0 
BETA (I) =SQRT(M(I) **2.0-1.0) 
CONTINUE 
CALL SHOCK({Z,M, THETA) 
CALL MAX (Z,M, DELTA) 
CALL AVG(Z,M,NU) 
CALL DISYSB(Z, DELTA, YSB) 
CALL CENTER (Z, THETA, BETA, SIGMA) 
CALL AREA(Z,M,A) 
CALL PRESS (Z,SIGMA,M, P) 
CALL SPRATIO(Z,P,A,NU, RATIO) 
CALL DISXO(Z, YSB, BETA, RATIO, THETA, X0) 
The following do loop, for each Mach number, calculates 
the angle of the asymptote of the hyperbola that 
defines the shock. Then the point on the x and y 
axis where the difference between the shock wave 
angle and the asymptote angle is less than 1 per 
cent is found. This defines the limit of the shock wave 
for purposes of calculating the entropy change across 
the shock. The change in entropy is then calculated. 
DO 15 J=1,2Z 
B(J) =X0(J) /BETA(dJ) 
CANGLE (J) =ATAN (B(J) /X0 (J) ) 
CALL SHAPE (X0(J),BETA(J) ,X,Y,) 
CALL SWANGLE (X0 (J) ,BETA(dJ) , ¥, SANGLE, ANGLED) 
DO 10 I=2,7500 
IF (SANGLE(I) .EQ.0) GO TO 10 
DIFF=(SANGLE (I) -CANGLE (J) ) /CANGLE (J) 
IF(DIFF.GT.0.01) GO TO 10 
IF(DIFF.LE.0.01) THEN 
CALL TOTAL(M({J) ,X0(J) ,BETA(J) ,X(I) , TOTALS) 
GO TO 15 
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END IF 


10 CONTINUE 
15 CONTINUE 


+ 
+ 
+ 


STOP 
END 


Calculate the shock wave angle for sonic flow behind 
the shock 
INTEGER I,2Z 
REAL M(Z), THETA(Z), GAMMA 
GAMMA=1.4 
DO 05 I=1,2 
THETA (I) =ASIN(SQRT((1.0/4.0*GAMMA*M(I) **2.0)) 
* ( (GAMMA+1.0) *M(I) **2.0-(3.0- GAMMA) 
+ (SQRT ( (GAMMA+1.0) * ( (GAMMA+1.0) *M{I) **4.0-2.0 
* (3.0-GAMMA) *M(I) **2.0+ (GAMMA+9.0))))))) 


05 CONTINUE 


05 


05 


RETURN 
END 


Calculate the maximum deflection angle, in radians, 
for a given freestream Mach number in air. For this 
calculation, the maximum deflection angle without shock 
detachment was calculated from NACA report 1135, chart 
5. An equation for the line of maximum deflection was 
obtained by performing a curve fit of the data 
presented in the chart. 
SUBROUTINE MAX(Z,M, DELTA) 
NLOG = Natural log of Mach number between 1.1 and 10 
TVAR = Natural log of shock wave detachment angle in 
degrees 
INTEGER I,2Z 
REAL M(Z), NLOG(10), TVAR(10), DELTA(Z) 
DO 05 I=1,2 
NLOG (I) = LOG(M(T) ) 
TVAR (I) =3.80568+0.41631*LOG (NLOG(T) ) 
DELTA (I) =(EXP(TVAR(I)))/57.29578 
CONTINUE 


Calculate the angle between the sonic line and a line 
normal to the freestream direction 
SUBROUTINE AVG(Z,M,NU) 
INTEGER I,Z 
REAL M(Z), NU(Z) 
DO 05 I2#1,2 
NU (I) =(3.9901+56.9298*LOG (M(Z))) /57.29578 
CONTINUE 
RETURN 
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END 


Calculate the distance from the body centerline to the 

sonic point on the body 

SUBROUTINE DISYSB(Z, DELTA, YSB) 

DIST1i=Distance along a line drawn tangent to the body 
sonic point from the x-axis to the sonic point on 
the body 

R=Radius of the body apex (in feet) 

ANG3=The angle between the x-axis and a line 
perpendicular to DIST1 at the body sonic point and 
extending from the sonic point to the x-axis 

INTEGER I,2Z 

REAL PI, YSB(Z), DIST1({10), R, ANG3(10), DELTA(Z) 

PT=3.14159 

_ R=1.5 

DO 05 I=1,2 

ANG3 (I) =PI- (PI/2.0) -DELTA(I) 
DIST1 (I) =R* ( (SIN(ANG(I) ) /SIN(DELTA(I) ))) 
YSB(I) =DIST1 (I) *SIN(DELTA(T) ) 
05 CONTINUE 
RETURN 
END 


Calculate the angle of the shock wave at the center 

of mass streamline between the shock wave vertex 

and the sonic point of the shock wave 

SUBROUTINE CENTER (Z, THETA, BETA, SIGMA) 

INTEGER I, Z 

REAL THETA(Z), BETA(Z), SIGMA(Z) 

DO 05 I=1,2 

SIGMA (I) =ATAN({ ( (TAN (THETA(I) ) ) **2.0* (3.0/2.0) **2.0 

+ - (1.0/BETA(TI) ) **2.0)*((3.0/2.0) **#2.0-1.0)) **0.5) 
05 CONTINUE 

RETURN 

END 


Calculate the area ratio for isentropic flow 
SUBROUTINE AREA(Z,M,A) 
INTEGER I,Z 
REAL GAMMA, M(Z), A(Z) 
GAMMA=1.4 
DO 0S I=1,Z 
A(I) =((GAMMA+1.0) /2.0) ** ( (GAMMA+1.0)/(2.0 


+ *GAMMA-1.0))) *M(I) *(1.0+( (GAMMA-1.0) /2.0) 

+ *M(I) **2.0) ** (- ( (GAMMA+1.0) /(2.0* (GAMMA-1.0)))) 
05 CONTINUE 

RETURN 

END 
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Calculate the ratio between stagnation pressure before 
and after the shock wave for a given shock wave angle 
SUBROUTINE PRESS (Z,SIGMA,M >) 
INTEGER I, Z 
REAL GAMMA, SIGMA(Z), M(Z), P(2) 
GAMMA = 1.4 
DO 05 I=1,2 
PI (I) =1.0/(( ( (GAMMA+1.0) *M(I) **2.0* 
SIN(SIGMA(I)) **2.0)/( (GAMMA-1.0) * 
M(I) **2.0*SIN(SIGMA(I))**2.0+2.0) ) ** (GAMMA/ 
(GAMMA-1.0) ) * ( (GAMMA+1.0) /(2.0*GAMMA*M (I) **2.0* 
SIN (SIGMA(TI) ) **2.0- (GAMMA-1.0)))** 
(1.0/ (GAMMA-1.0))) 


++ + + + 


05 CONTINUE 


05 


05 


RETURN 
END 


Calculate the ratio of the distance from the x-axis to 
the shock sonic point to the distance from the x- 
axis to the body sonic point 
SUBROUTINE SPRATIO(Z,P,A,NU,RATIO) 
INTEGER I,Z 
REAL P(Z), A(Z), NU(Z), RATIO 
DO 05 I=1,2 
RATIO(I)= 1.0/((1.0-P(I) *A(I) *COS (NU(I))) **0.5) 
CONTINUE 
RETURN 
END 


Calculate the distance from the foremost point of the 
detached shock wave to the intercept of its asymptote 
on the x-axis 

SUBROUTINE DISX0(Z, YSB, BETA, RATIO, THETA, X0) 

INTEGER I, 2Z@ 

REAL YSB(Z), BETA(Z), RATIO(Z), THETA(Z), X0(Z) 

DO 05 I=1,2 

X0 (I) =YSB(I) *BETA(I) *RATIO(TI) * ((BETA(I) ) **2.9* 

+ (TAN (THETA(I) )) **2.0-1.0) **0.5 

CONTINUE 

RETURN 

END 


Calculate the x and y coordinates for any point on the 
shock wave 
SUBROUTINE SHAPE (X0, BETA, X, Y) 
INTEGER I 
REAL X(7500), Y¥(7500), BETA, XO 
X(1)=0.1 
Y¥(1)=0 
DO 05 T=#2,7500 
X(I) =X(I-1)+0.2 
IF(X(I).LT.X0) THEN 
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Y¥ (I) =0 

ELSEIF (X(I).GE.X0) THEN 

Y (I) =SQRT (X(I) **2.0-X0**2.0) /BETA 
END IF 

05 CONTINUE 

RETURN 


END 


Calculate the angle between the freestream and a line 
drawn tangent to the shock wave at any point 
SUBROUTINE SWANGLE (X0, BETA, Y, SANGLE, ANGLED) 
INTEGER I 
REAL X00, BETA, Y(7500), SANGLE(7500), PI, ANGLED(7500) 
PI=3.14159 
DO 05 I=1,7500 
IF (Y(I) .EQ.0) THEN 
SANGLE (I) =0 
ELSEIF (Y(I).GT.0) THEN 
SANGLE (I) =ATAN ( (SQRT (X0**2.0+BETA**2.C*Y(I) **2.0) ) 
/ (BETA**2.0*Y(TI))) 
ANGLED (I) =SANGLE (I) *180/PI 
END IF 


05 CONTINUE 


RETURN 
END 


Calculate the total change in entropy across the shock 

wave. All calculations are in SI units 

RO=The density of air. For these calculations this 
value is taken from the ICAO standard atmosphere 
tables at an altitude of 30,000 meters. 

TEMP=The temperature in degrees Kelvin of the 
atmosphere at 30,000 meters altitude, taken from 
the ICAO standard atmosphere tables 

GAMMA=Ratio of specific heats 

SGC=Specific gas constant. Assuming an ideal gas, for 
air this value is 287 Joules/Kg-K 

SSA=Local speed of sound in air 

VEL=Velocity, in meters per second, of the freestream 

SHCV=Specific heat of air at constant volume 

PI=Radian equivalent of 180 degrees 

LLINT=Lower limit of integration. The starting point, 
located at the shock wave vertex, for numerical 
integration to calculate the entropy change across 
the shock wave 

ULINT=Upper limit of integration. The ending point for 
numerical integration across the shock wave 

H=The integration step size for application cf 
Simpson’s rule 

N=The number of Simpson’s rule integrations, and the 
ending point for summation of individual 
integrations 
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CENTERX=The center point, along the x-axis, for each 
Simpson’s rule integration 


TSTEP=Specifies the beginning (TSTEP(1)), center 
(TSTEP(2)), and end (TSTEP(3)) points for 
Simpson’s rule application for each individual 
summation 


RADIUS=The distance from the centerline x-axis to the 
shock wave at each TSTEP location 
TANGLE=The shock wave angle at each point of 
intersection of the radius and the shock wave 
DELTAE=Change in entropy per mass flow across the shock 
at each point of intersection of RADIUS with the 
shock wave 
DELTAS=Distance along the shock at each interval of 
integration 
HOLD=Temporary storage variable 
INTEGER I, K 
REAL XO, X, BETA, DELTAE(3), TOTALS, LLINT, ULINT, H, 
+ CENTERX, PI, RADIUS(3), SHCV, GAMMA, TANGLE(3), 
+ TSTEP(3), RO, VEL, SSA, TEMP, M, DELTAS 
RO=.01841 
TEMP=226.509 
GAMMA=1.4 
SGC=287 
SSA=SQRT (GAMMA*SGC* TEMP) 
VEL=M*SSA 
SHCV=717.5 
PI=3.14159 
LLINT=X0*.3048 
ULINT=X* . 3048 
HOLD=0 
H=0.5 
N=INT( (ULINT-LLINT) / (2.0*H) ) 
CENTERX=LLINT+.05 
DO 0S I=LLINT, N 
DO 10 K=1,3 
TSTEP (1) =CENTERX- .05 
TSTEP (2) =CENTERX 
TSTEP (3) =CENTERX+.05 
IF (TSTEP(K) .LE.LLINT) THEN 


RADIUS (K) =0 
TANGLE (K) =0 
DELTAE (K) =0 
GO TO 10 
END IF 


RADIUS (K) = (SQRT (TSTEP (K) **2.0-LLINT**2.0))/BETA 
TANGLE (K) =ATAN ( (SQRT (LLINT**2 .0+BETA**2 .0*RADIUS (K) 


+ **2.0))/ (BETA**2.0*RADIUS (K) ) ) 

DELTAE (K) =SHCV* (LOG ( (2. 0*GAMMA*M* *2 .0*SIN (TANGLE (K) ) 
+ **2.0- (GAMMA-1.0)) / (GAMMA+1.0) ) -GAMMA* 

+ LOG ( ( (GAMMA+1.0) *M**2.0*SIN (TANGLE (K) ) **2.0) / 
+ ( (GAMMA-1.0) *M**2.0*SIN (TANGLE (K) ) **2.0+2.0))) 
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10 CONTINUE 
DELTAS= (0.1*.3048) /COS (TANGLE (2) ) 
HOLD=HOLD+ ( ( (H/3.0) * (DELTAE (1) +4.0*DELTAE (2) 
+ +DELTAE (3) )) *RO*2.0*PI*RADIUS (2) *DELTAS*VEL* 
+ SIN (TANGLE (2) ) ) 
CENTERX=CENTERX+0.1 
05 CONTINUE 
TOTALS =HOLD* TEMP 
RETURN 
END 
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